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ABSTRACT 

We describe an analytic distribution function of a finite, oblate stellar system that is 
useful for the practical modelling of dark halos. The function is determined by lowering 
Evans's (1993) distribution function of a fiattened, cored isothermal system in analogy 
to the lowering of the singular, isothermal sphere in the definition of the King (1966) 
models. We derive analytic expressions of the density, maximal streaming velocity and 
velocity dispersion profiles as a function of the potential. As for the King models, 
the potential must be calculated numerically. We also present a recipe for generating 
N-body realizations of this distribution function and examine the stability in three 
models with dimensionless spins A =0.0, 0.05 and 0.18 using N-body simulations with 
50,000 particles. The A = 0.18 model is unstable to the formation of a triaxial bar 
within ~ 5 King radii while the other models appear stable. We conclude that the 
slowly rotating systems are useful for modelling fiattened dark halos. 

Subject headings: galaxies: kinematics and dynamics - galaxies: structure - galaxies: 
Elliptical - methods: numerical 

1 Introduction 

The dark halos surrounding disk galaxies are probably not spherical. In order to simulate 
the effects of non-spherical halos on disks, it is useful to have an analytic distribution function 
(df) for flattened halos. Such models do exist, but they are either not very realistic models 
for galaxy halos (such as the flattened Plummer models of Lynden-Bell 1962) or inflnite 
in extent (for example, the scale-free logarithmic models of Toomre 1982). In this paper 
we construct a set of analytic DF's for flnite-radius, oblate galactic halo models, and show 
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that they can be sampled efficiently for use as an equilibrium starting condition in N-body 
simulations. 

Distribution functions are most useful when specified as functions of integrals of motion. 
Most realistic galaxy potentials admit three of these, but in the case of axisymmetric models 
explicit expressions only exist for two (except for the special cases of Stackel potentials and 
of spherical models). The DFs that we construct here will therefore be functions only of 
the two classical integrals of motion, namely the particle energy E and angular momentum 
about the symmetry axis L^. As shown by Lynden-Bell (1962), there is a unique relation 
between the even part of the DF, /(£", |iv^|), and the density p{R^ ^) in the meridional plane 
when expressed as a function of the potential and cylindrical radius. Arbitrary terms odd 
in Lz may be added to the DF without affecting the density p: since such terms effectively 
change the azimuthal direction in which some orbits are traversed, they set the total angular 
momentum of the system. 

In this paper, we modify an existing analytic DF for infinite fiattened halo models, due 
to Evans (1993), making it finite in extent. The resulting models have a range of fiattenings, 
central densities, outer radii and concentration parameters. Evans's model is described in 
§2, and our modification of it in §3. Section 4 contains a recipe for generating N-body 
realizations of the models, and the results of N-body experiments to test the stability of 
several of the models. We find that as long as the mean rotation of the models is kept 
within reasonable cosmological expectations, there is no sign of large-scale bar instability. A 
summary is provided in §5. 



Evans (1993) found the exact two-integral distribution function (df) for all axisymmetric 
Binney (1981) potentials, by applying Lynden-Bell's (1962) Laplace transform method. His 
models are among the few fully analytic axisymmetric potential-density-DF sets that can be 
written down, and they are very simple, at that: Binney 's potential is 



2 Evans's distribution function for the Binney potentials 





the mass density is 



p{R,z) 
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and the two-integral DF has the form 



f{E, L,) = All exp(-2^/a2) + B exp(-2^/a2) + Cexp(-^/a2). 



(3) 



The constant i?i (which Evans set to one) is introduced here in the potential for consistency 
of units. The related density. 



LOWERED EVANS MODELS 



3 



serves as a density scale parameter for the models. The constants A, B and C depend on 
the velocity dispersion ctq, the axis ratio of the equipotentials q and the core radius i?c, 
as follows: 

, 8(1 - q')Gpl mlGpl i2q' - l)p. 

In particular, spherical models have A = and coreless models have B = 0. The familiar 
isothermal sphere DF is recovered when both these constants are zero: in that case, the DF 
is a Maxwellian with velocity dispersion ctq and density pi exp( — ^/cTq). The coreless, B = 
models were first discovered by Toomre (1982). 



3 Lowered Evans models 



All the models (3) are infinite in extent; moreover, at large radii the mass density falls 
as inverse-square radius, implying that they have infinite mass. These models are useful 
for the purposes of modelling galactic halos, since they have an adjustable core radius and 
a fiat rotation curve at large radii. Nevertheless, for many applications it would be more 
convenient to have a finite-mass model, which nonetheless has a fiat rotation curve over an 
appreciable radius range. In analogy with the lowered isothermal, or King (1966) models, 
we were therefore motivated to construct finite oblate halo models by 'lowering' Evans's DF, 
in effect imposing a maximum energy on the stars in the model: thus we write 

r [{ALl + B)exp{-E/a'o) + C][exp{- E / a'^) - 1] if ^ < 0, 
fiE,L,) = l (6) 

V otherwise. 

Only stars with negative energy are included in these models. There are many other ways 
of accomplishing an energy cutoff in the DF, or indeed of limiting the DF in a different 
way (Binney and Tremaine 1987, Kashlinsky 1988, Rowley 1988) but we follow King's lead. 
The outermost radius at which the potential is negative (and hence the density non-zero) is 
usually called the tidal radius. 

The lowered DF (eq. 6) no longer corresponds to a self-gravitating system in Binney's 
potential (eq. 1). To contract self-consistent models with this DF, we use the fact that the 
DF determines the mass density in terms of the gravitational potential: 

dvRdv^dv, f = — , dLJE f = /j(*,i?). (7) 

Therefore we can achieve self-consistency by combining this expression with the Poisson 
equation, to yield 

V'* = 47rG'/j(*,i?), (8) 

and solving for the potential. 
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In the case of the lowered Evans DF, some algebra shows the corresponding density in 
an arbitrary potential (see eq. 7) to be 



= ]^TT'^l^al{AR^al + 2B) erf (V-2*/(7o) exp(-2*/(7, 



+ (27r)3/V3(C -B- AR^al) erf (V^/ao) exp(-*/a2) 
+7rV-2*[(7^(3A(7^i?2 + 25 - 4C) + ^*(2C - Aa^i?')], (9) 

where erf(a;) = 27r~^/^ Jq^ exp(— t^) is the usual error function. The maximum streaming 
velocity is obtained by inserting an extra factor of \Lz\l R in the integral (7): then 

/ot^^max = 27r*2(Ai?V2-C) + 27r*(72(2C-5-3Ai?V2) 

+TTal{lAR^al + 3B - 4C) + iira^iC -B- 2Ai?V^) exp{-E/al) 
+Tr(j^{AR^(j^ + B)exp{-2E/(j^). (10) 

The velocity dispersions follow similarly. They are given by af = J d^vfvf: 



^7rij^^f^[AR^a^o - 2C] + 27rija^o^f^[-AR^a^o " + ^ 



1-7 

+Tr4^l^[-AR'a'o + W- 4C] 



+ (27r)3/Vo^[-Ai?V2 -B + C] exfi^/ao) exp{-tP/a',) 

+ i7r3/Vo'[Ai?V2 + 2B] edi^f^/ao) exp(-2V'/a2), (11) 



anc 



pal = ^TrV-^y^iaAi^Vo^ - 2C] + 27r,l;a',,/^^[-3AR'a', -h + ^C] 

.21 



+Tra*^^-2ij[—AR^al + 3B - 4C] 
+ (27r)3/Vo^[-3Ai?V2 - 5 + C] erf(y^/ao) exp{-tP/a',) 
+ ^7r3/Vo'[3Ai?V2 + 25] erf(y^/ao) exp(-2V'/a2). (12) 

At i? = all three dispersions are the same (as they are in all non-singular two-integral 
models). The difference pa^ — paj^ indicates the anisotropy of the velocity dispersion. This 
expression is proportional to AR^cTq (the case A = corresponds to spherical models, since 
the B- and C-terms are isotropic), with a coefficient that is positive for all ^ < 0: therefore 
all oblate non-rotating models (A > in eq. 3) satisfy > (Tr everywhere, while the 
converse holds for the prolate models. 

The King (1966) models are just the special case A = i? = 0. These models are spherical, 
with Rc = 0. (Note the distinction between the core radius Rc of the Binney potential before 
lowering the DF, and the 'King' radius of the density, which measures the core radius of the 
self-consistent mass density of the model. Both these radii lose their meaning somewhat in 
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the lowered models.) In the spherical case it is possible to solve Poisson's equation (eq. 8) 
as a simple forward integration in radius. Given choices for the density scale pi and the 
velocity dispersion parameter ctq, the central potential energy (we have fixed the cutoff 
energy at zero) may be specified as a boundary condition for this integration. Deeper central 
potential depth implies a higher central density and a larger tidal radius, and therefore leads 
to a more centrally condensed model. 

Calculating the two-dimensional models is a little harder, since we cannot now integrate 
the Poisson equation directly. However, the principle remains the same: once we have chosen 
the scale parameters ctq and Rc and the fiattening every value of the central potential (or 
density) implies a unique non-singular self-consistent model. In practice, numerical iterations 
are required to obtain the density and potential: a guess is made at the potential, the density 
corresponding to the DF in this potential is calculated from eq. 6, Poisson's equation is solved 
for the potential corresponding to this density, and this new potential is taken as the start of 
the next iteration. In this work we have used a multipole expansion (Prendergast & Tomer 
1970) to solve for the potential: the models are quite smooth, and hence it was profitable 
to solve for the radial dependencies of the dominant harmonic terms (a few one-dimensional 
functions) rather than attacking a two-dimensional grid calculation. Eor finite-extent systems 
the multipole expansion has the further advantage that boundary conditions at infinity are 
automatically handled correctly. To ease numerical convergence, the higher harmonics were 
introduced one at a time, allowing a few iterations for any oscillations to stabilize before the 
next term was added. Potentials of the models were calculated up to / = 4, and the densities 
were then derived using eq. (9). 

King models have non-singular cores, in spite of the fact that the DF does not contain 
the i?-term (the only term proportional to the Binney potential core radius Rc). Eor very 
deep central potentials these models have sharply peaked densities, though, asymptoting 
to the singular isothermal sphere in their central regions. The central core radius, or King 
radius, of the King models' potential is related to the central density pc by 

At this radius the gravitational potential has risen by about 2(Tq over the central value, 
provided the potential well depth is well below ctq. A concentration parameter is usually 
defined as the ratio of the outer radius of the model and the King radius. 

The fiattened models with B = also become more centrally concentrated as the depth 
of the potential is increased. However, these models are not very satisfactory: the central 
density contours have undesirable depressions on the symmetry axis. It turns out that a 
small i?-term (i.e., non-zero Rc) suffices to remove this effect, since the addition of this 
isotropic term, more centrally concentrated than the C-term, serves to round off the central 
density figure. Some experimentation shows that the choice Rc = vk results in a more 
elliptical central core, but does not affect the fiattening of the halo model at larger energies 
too drastically. Addition of a i?-term affects the central density: in practice, a quadratic 
equation in R^ needs to be solved to find the core radius that yields a central density whose 
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Fig. l.-Iso- density contours for two q = 0.8, $o = — ScTq models. The density on adjacent contours 
differs by a factor of two. The left-hand panels show the model with i?-term (which has King radius 
0.053), the right-hand panels the one without (in this case = 0.163). Both models have similar 
shape outside the core, but the addition of the i?-term removes the strong 'peanut' dimples in the 
central regions. 
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King radius is equal to i?c- The isodensity contours for q = 0.8, = ~6(Tq models with 
and without such a i?-term are shown in Eigure 1, illustrating the rounding effect in the 
central regions. Eor shallower models (^o~ — 3(Tq), it is better to take Rc = krx with A; < f , 
otherwise all the models end up rather round. 

Our models, in summary, are the self-consistent realizations of the axisymmetric distri- 
bution functions (6). The parameters ctq, /Ji, and q may be chosen freely, after which the 
constant Rc is adjusted until the central density pc is equal to {*^k^ (Tq j ^tt G R^^) ^ i-^- until the 
King radius is equal to Rc/k^ where k lies between 0.3 and f . 

The choice of the four parameters ctq, q and is equivalent to picking a tidal 
radius, a concentration parameter, a flattening and a central density (or a total mass) for 
the system. Eigures 2 and 3 present some of the relevant relations: they show the dependence 
of the concentration parameter rt/rx and of the scaled central density as a function of the 
scaled potential well depth, ^q/cq- 

Eigure 4 shows the circular velocity curve in the equatorial plane of the models of 
Eigure f. 

4 N-BODY Realizations 

The lowered Evans model is useful for examining galactic dynamical problems that depend on 
flattened potentials such as disk warping or polar rings. The model is finite in extent making 
it ideal for use in N-body simulations. In this section, we present a recipe for generating 
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Fig. 4. -Circular velocity curves in the equatorial plane of the q = 0.8, $o = ~6(Tq models. 
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an N-body realization of a lowered Evans model. We then test the stability of three sample 
models using q = 0.8 with different amounts of spin: a non-rotating model, a maximally 
streaming model and a model with spin corresponding to the cosmological expectations of 
the dimensionless spin parameter, A = 0.05. 

4.1 Initial Conditions 

We can generate an N-body realization from any DF of the form f[E^Lz) by sampling 
from it in two stages. Eirst, we sample values of R and z from the density distribution, 

z) to find the particle positions. Then, for each position the DF and the gravitational 
potential define the distribution of velocities, and we sample from this function to assign 
each particle's velocity. 

We use the acceptance-rejection technique for sampling the distributions (e.g.. Press 
et al. f993), which works as follows. Let -Fmax be the maximum value of the distribution 
function. Then, we randomly select, with a uniform distribution, a point in the allowed 
domain of the independent variable(s), and we also sample a 'test' value F from the range 
[0, -Finax] of the distribution function. The value of F determines whether the sampled point 
will be included or not: if F is less than the value of the distribution function at the sampled 
point, that point is accepted as valid, otherwise it is rejected and another point is sampled. 
If the distribution being sampled is very non-uniform, this process may be quite inefficient, 
with many of the sampled points being rejected. To avoid this situation, it is worthwhile 
to transform the independent variables to coordinates whose distribution is as uniform as 
possible. 

The distribution of the particle positions (up to a normalization) is the density p{R^ z). 
Eor the lowered Evans model, p ~ r~^, so it is convenient to introduce the variable u = 
tan~^ zj R which makes the mass element almost uniform: 

p{R, z)dRdz = p{u, R){R^ + z^)dRdu. (14) 

The domains of R and u are < i? < rj and — 7r/2 < u < Tr/2. Eor every (u^R) point 
sampled, we also sample an independent random azimuthal angle ^, finally allowing us to 
define x = Rcos cf), y = i?sin cf) and z = i?tan u. Each particle thus sampled is then assigned 
a velocity vector, in essentially the same way: the velocity distribution function at each 
particle's position is known from the DF (eq. 6), where F = -|- ^(i?, z) and L = Rv^, and 
the velocity is always less than the local escape speed, Vesc = a/— 2^. The velocity vector can 
thus be sampled with the acceptance-rejection technique from inside a sphere of this radius. 

The DF of eq. (6) does not depend on the sign of L^: consequently it has no net 
streaming. We can introduce azimuthal streaming into the model by varying the number of 
particles with positive and negative. In a non-rotating model, there are equal numbers 
of particles with going in opposite directions while for a maximally streaming model, the 
sign of is the same for all the particles. All intermediate cases have varying fractions of 
going in opposite directions. 
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4.2 Stability 

Although it is easy to formulate an equilibrium DF, it is not guaranteed that the resulting 
system is stable. We therefore tested the stability of rotating and non-rotating models with 
50,000-particle N-body simulations. Initial conditions were generated with the procedure 
described above (Figure 5). We worked in units in which G = 1 and pi = (47r)~^, and 
investigated models with ctq = 2~^/^, flattening q = 0.8 and central potential = — 6.0(Tq. 
The parameter Rc was chosen to be equal to the King radius, as described in §3. The base 
model has the following properties: (i) The central density is pc = 126 so that the equivalent 
King radius from equation (13) is vk = 0.053, (ii) the model extends to a "tidal" radius 
rt = 2.14, (iii) the core crossing time is Tcore = (Stt/G/Jc)^''^ = 0.30, and (iv) the system 
crossing time is Tsys = 2rt/(To = 6.0. We introduce varying degrees of rotation in three 
models by varying the fraction of particles going in opposite directions. We parameterize 
the rotation using the dimensionless spin parameter, A = L\E\^^'^ M~^^'^ ^ that is used 
to quantify spin in dark halos formed in cosmological models. The typical value for a 
cosmological dark halo is thought to be A = 0.05 (e.g. Barnes & Efstathiou 1987; Warren et 
al. 1992). The three models we investigated have A = 0.0 (a non-rotating halo), A = 0.05 (a 
"cosmological" halo) and A = 0.18 (the maximally streaming halo). 

We then used a tree code (Barnes & Hut 1986; Dubinski 1988) to simulate the models 
for 24 units of time corresponding to 80 core crossing times and 4 system crossing times. 
We used an opening angle tolerance, ^ = 1.0 and calculate cell-particle forces to quadrupole 
order with a particle softening radius, r^o/t = 0.005. We integrate the trajectories using a 
leapfrog integrator with a timestep At = 0.02 corresponding to 12 steps per core crossing 
time. The error in the total energy of the systems was no more than 2.5% by then end of 
each run. 

We measured the stability of the models by comparing the density and velocity disper- 
sion profiles, averaged over spherical shells, at early and late times. Shell averages are trivial 
to compute from N-body simulations, requiring a simple binning of the particles by radius. 
The exact shell averages that correspond to the initial analytic DF are also straightforward 
to obtain: one can show that for a given cylindrically symmetric function g{R^ z) the spher- 
ically averaged function 'g{r) found by averaging the function g within a thin spherical shell 
at radius r is 

g{r) = - r g{[r'-z']'/\z)dz. (15) 
r Jo 

We can therefore calculate the averaged density profile 'p{r) and velocity dispersion profiles 
(t|j, cr^, and (t| by inserting the functions from equations (9) and (11), 12) into this integral. 

Figure 6 compares the expected averaged density profile to measurements from the 
rotating and non-rotating models at the final time t = 24.0. The density profiles remain 
unchanged over a large range in radii except near the center of the models. The density 
declines slightly in the center, probably in response to the greater degree of integration error 
within the core. Nevertheless, the mass profiles of the models appear stable over at least 
four system crossing times. 
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t= 0^0 



Fig. 5. -Edge on view of the 50,000 particle model used in the simulations. The box width is 4.0 
units across (80 King radii). The flattening in the density is about qp = 0.6 in the center increasing 
to qp = 0.8 at the tidal radius. 
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Fig. 6. -The spherically averaged density profiles for the 3 models normalized to the central density. 
The profiles of the different models are offset by one dex for illustration. The curves are the 
theoretical expectation and the points are the density at the end of the simulation, estimated by 
binning particles in spherical shells. The agreement is very good over a large range in radii except 
for a slight dip in the core that is probably due to integration error. 
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Fig. 7a. -Axial ratio profiles of density contours at t = 0.0 and t = 24.0 for the three models. The 
profiles for the A = 0.0 and A = 0.05 models remain fairly constant. There is a hint of a bar forming 
within r = tk though it is difficult to judge if this is real because of integration and measurement 
error. The A = 0.18 model shows a distinct bar out to r ~ Sr^ indicating the onset of the bar 
instability. 

A more refined indication of stability is provided by a study of quadrupole terms in 
the mass distribution. We therefore estimate the axial ratio profile of the dark halos using 
the technique described in Dubinski & Carlberg (1991). In this algorithm, initial values for 
the axis ratios qi and q2 are assumed, and used to calculate a starting approximation to the 
modified inertia tensor, li^ = XiXj/a'^ for particles in ellipsoidal shells of axial ratios qi and 
q2 {xi is the particle position and = + / qf + z"^ I q^ is the particle elliptical radius). 
New axial ratios, and the orientation of the ellipsoidal figure, are then estimated from li^ 
through q\ = lyyjlxx and q^ = Izz/Ixx^ and used to calculate an improved approximation 
to the modified inertia tensor. Starting with particles in a spherical shell [qi = q2 = 1), 
this process is repeated for several iterations until the axial ratios and the shell orientation 
converge to values within a specified tolerance [Aq = .001). 

Eigure 7 presents the axial ratio profiles measured from the simulations at the initial 
and final times. The isodensity contours of the dark halo are slightly peanut shaped near 
the center so that the estimate of the axial ratio from the modified inertia tensor (which 
assumes that the density contours are ellipsoidal) will underestimate the ratio of the extent 
of the isodensity contours along the R and z axes somewhat. The axial ratio profile does not 
change dramatically for the A = 0.0 and A = 0.05 models suggesting that these systems are 
not susceptible to strong bar instabilities, though a careful look at the A = 0.0 and A = 0.05 
halos reveal a small hint of a bar within one King radius. Since it did not extend far beyond 
the core, it is difficult to say whether this feature is real, or a result of integration error 
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t= 



t = 24 



Fig. 8. -Particle plot of the center of the A = 0.18 as viewed down the 2;-axis at the t = 0.0 and 
t = 24.0. The length of the box is / = 20rfc. At t = 24.0 a triaxial bar is visible in the center of the 
system. 

and particle noise. On the other hand, the A = 0.18 model is clearly unstable, forming a 
distinct triaxial central bar with an axis ratio qi ~ 0.7 and q2 — 0.5 within 5 King radii. 
A particle plot at the end of the simulation (Eigure 8) show this bar clearly: inspection 
of earlier snapshots reveals that the bar grows during the first system crossing time. It is 
not too surprising that the A = 0.18 model is bar-unstable, since this type of instability 
generically affects fast-rotating systems: though it is perhaps remarkable that even a system 
as hot as this one cannot quench it. 

The spherically-averaged velocity dispersion profiles (Eigure 9) provide further evidence 
for the stability of the slowly-rotating models, though there is a slight deviation within a 
King radius, again probably refiecting some integration error within the core. The instability 
in the A = 0.18 model, in contrast, shows up clearly as it develops a noticeable dip in the cr^ 
profile. 

In conclusion, the lowered Evans models appear to be stable when the rotation rate is 
small, but they may suffer bar instabilities in the extreme case of maximal streaming. We 
recommend the non-rotating and cosmological halos for application to halo modelling but 
urge caution in the use of more rapidly rotating models. In any case, the bar instability 
becomes apparent within a crossing time and can be checked in practice before applying the 
model to a problem. 
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Fig. 9a. -Spherically averaged velocity dispersion profiles for the three models. The curves rep- 
resent the theoretical shell averaged profiles while the points are the dispersion of the R, cj) and z 
components of velocity at the end of the simulation, measured by binning particles in shells. The 
profiles remain unchanged for the A = 0.0 and A = 0.05 models over most radii though there is 
a slight deviation within the core again a refiection of integration and measurement error. The 
systematic dip in in the A = 0.18 model is another manifestation of the bar instability. 
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Fig. 9c.- 
5 Summary 

The lowered Evans distribution function provides equilibrium oblate stellar systems 
which are a flattened analogue of the spherical King models. The density and velocity dis- 
persion profiles can be expressed analytically in terms of the potential though the calculation 
of the potential still requires a fair amount of numerical computational effort. The payoff 
for this effort is a finite model ideal for N-body simulations of galactic systems involving 
fiattened dark halos. Eurthermore, we provide a simple recipe for generating an N-body 
realization of the distribution function. A sample of N-body simulations shows that the 
models are stable for slowly rotating models with spin corresponding to cosmological dark 
halos. The maximally streaming model is unstable to bar formation despite its high dy- 
namical temperature. We therefore caution users of these models to watch out for the bar 
instability in rapidly rotating models. In the near future, we plan to apply these models to 
the formation of warps in disk/halo systems. 
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